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Abstract 


The derivation of gravity anomalies at the surface of the earth from 
satellite to satellite tracking or satellite gradiometry observations is a down- 
ward continuation problem. In developing the spectral interrelations in terms 
of spherical harmonics between different gravity quantities such as the disturb- 
ing potential T, gravity anomalies Ag, and i_| at the surface and in satellite 

Sr 2 

altitude a discussion of the features of these two methods for the estimation of 
surface gravity information is possible. 

For a local improvement of our knowledge about the gravity field inte- 
gral formulas in the parameter domain have to be applied instead of a 
representation by spherical harmonics. The neglected regions will cause a 
truncation error. The application of the discrete form of the integral equations 
connecting the satellite observations with surface gravity anomalies is discussed 
in comparison with the least squares prediction method. 

One critical point of downward continuation is the proper choice of the 
boundary surface. Practical feasibilities are in conflict with theoretical 
considerations. The properties of different approaches for this question are 
analyzed. 

As a result the considerations indicate the possibility of deriving mean 
gravity anomalies at the surface of the earth. By taking into account theoretical 
restrictions these anomalies are comparable with terrestrial gravity anomalies. 
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Introduction 


By the "Earth and Ocean Physics Applications Program"(1972) from NASA 
the development of some new satellite techniques is supported that will enable the 
determination of the earth’s gravity field in a different way. Specifically, the 
methods of altimetry, satellite to satellite tracking, and satellite gradiometry are 
expected to lead to a considerable improvement in our current knowledge of the 
gravity field, hi several theoretical studies the capability of these techniques is 
analyzed and compared, for example by Kaula (1969), Glaser and Sherry (1971), 
and Forward (1971). Since the successful results of Muller and Sjogren (1968) for 
the recovery of the gravitational field of the front face of the moon by a direct 
mapping, satellite methods received growing attention not only for a global but 
also for a local analysis of the gravity field. Local solutions would considerably 
reduce the number of unknowns and can be concentrated in unsurveyed areas. 
Simulation studies by Reed (1973) for satellite gradiometry and by C.R. Schwarz, 
(1970) and Hajela (1974) for satellite to satellite tracking brought promising re- 
sults for this purpose. 

Chovitz (1973) turned new attention to the satellite analysis methods as a 
downward continuation problem. Especially for the local recovery serious trouble 
can be expected and the theoretical considerations for example by Moritz (1970), 
Bjerhammar (1973), and Krarup (1969) have great practical importance. 

This report has the intention to present a unified smoothing scheme by 
spectral analysis according to the idea of Meissl (1971). It shall allow us to com- 
pare different satellite techniques -he re satellite to satellite tracking and satellite 
gradiometry-by their sensitivity with respect to a certain unknown. The analysis 
of the downward continuation operator involved in the solution equations will 
illuminate the scale of difficulties connected with the evaluation of surface gravity 
anomalies due to this problem. Thereby the specific properties of the prediction 
and the deterministic approach will be compared. 

2. The spectral properties of the operators related to the earth's potential. 

For an optimal design of geodetic satellite experiments it is of basic im- 
portance to analyze the resolution and accuracy of the desired quantity based on 
the expected observation precision. The final aim will be to derive the equipoten- 
tial surface at mean sea level, the geoid. This surface may be computed from 
gradients, gravity disturbances, deflections of the vertical, gravity anomalies, 
or from the disturbing potential itself. The derivation of these quantities is 
possible by terrestrial measurements but also from different types of satellite 
data or from a combination of two or more of these techniques, hi order to 
find out the specific character of these different methods with respect to the 
desired quantity a unified model for their comparison would be of great use. 


1 


The intention of an unified model will not be to present a precise analysis 
for a special experiment but it shall serve a feeling of possibilities and limitations 
for the derivation of a gravity quantity from different types of observations. There- 
fore to a certain degree the conclusions will always suffer from simplifications 
introduced to obtain an easily understandable closed model. 

Because many properties of a physical quantity and its relations to others 
may be well described by spectral analysis we develop the desired unified model 
by expanding gravity quantities into spherical harmonics. An analysis based on 
this idea was initiated by Meissl (1971). 


In a first step we concentrate our interest on the interrelations between the 
disturbing potential T, the gravity anomalies Ag, and the first and second deriv- 
ative of the disturbing potential, and later on denoted as gravity quantities 

all of them on a sphere a with mean radius R of the earth. This collection is per- 
haps not representative but it has the advantage that these quantities are harmonic, 
their interrelations are simple, and their spectral properties reach from a high 
amplitude J,ow frequency spectrum for T to a low amplitude high frequency spec- 



The general form of a spectral representation by expansion into spherical 
harmonics for a gravity anomaly on ct is (Heiskanen-Moritz, p. 255) 


“A _ _ 

(2.1) Ag(P') = ) > (c(Ag) nn R na (P , ) + s(Ag) nn S nn (P / )) 

Z-j L-1 
n-0 m=0 

Ag(P # ) . . . gravity anomaly at a surface point P'ce 


c(Ag) no ) coefficients of the expansion of Ag into spherical harmonics of 
s(Ag) nn J * " ’ degree n and order m [in units of gravity] 

Rnn(P H . . . fully normalized spherical harmonics in p\ 

Snn (P ) I 


The relation between the disturbing potential T and the gravity anomalies Ag 
expressed by Stokes formula is 

(2. 2) T(Q^ = ^ J Ag(p') St((i))da 

a 

Expansion of the Stokes operator into spherical harmonics (Heiskanen-Moritz, 
p. 97 (2-169) and p. 33 (1-82')) leads to 
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oo a 

(2.3) St(p',Q') = £ -£TT £ (R n »(P / )R.g(^)+ S nB (P / )S nB (Q / )) 

n~ 2 m- 0 

Insertion of equation (2.1) and (2. 3) into (2. 2) and using the orthogonality relation- 
ships will give 

OO XI 

(2.4) T(QV£ (c(Ag) nB R nm (Q / )+s(Ag) nB S nB (Q / )), 

n= 2 m= 0 

a spherical harmonic expression for T(Q ; ). 

We may also develop T(Q X ) directly into spherical harmonics and get, 
similar as for Ag(P 7 ) in equ. (2. 1) 

00 n 

(2.5) T(Q> Y y’c(T) nm R ni! (Q / )+s(T) nB S nB (Q / )), 

— J — « 

n= 2 0 

with c(T) nn ) coefficients for the spherical harmonic representation of T 
s(T) nn / ' * * [in units of potential] . 

A comparison of the coefficients in equ. (2.4) and (2.5) shows 

(2.6) fc(T)n«Y R ( c(Ag) na 

\s(T)n „ J n_1 \s(Ag) n , J 

The same relationship holds true for the coefficients of an expansion into 
Legendre polynomials (Heiskanen-Moritz, p. 97) 

(2.7) T n = ~Y Ag n 

T n > Ag n . . c coefficients of the Legendre polynomial for T and Ag 
[units of potential and gravity] . 

Similar expressions can be derived for the other harmonic gravity quantities. 
From the spectral representation of their interrelating operators we find the least 
upper bounds (LUB) and greatest lower bounds (GLB) of the operators. For the 
quotient in equ. (2. 7) we derive for example 

TV 

GLB ( — =- ) = 0 for lim n -* 03 and 

LUB A) = R for n= 2. 
n-l' 
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The bounds are the result of the frequency properties of the operator, which in 
this example will damp high frequencies and reach its maximal amplitude in the 
low frequencies (maximum R for n= 2). The theoretical background of these 
considerations for different types of gravity quantities is described in Meissl 
(1971). 


The spectral interrelations of the quantities in mind and the GLB and 
LUB for the operators are collected in Table 1. 



T" 

! 

Ag„ 

(s). 

M 

T„ 


R/(n-l) 
0; R 

-R/(n+l) 
0; R/3 

R a /(n+l)(n+2) 
0; R s /12 

n 


1^91 


R(n-l)<n+l)(n+2) 
0; R/12 

©. 

' 

-(n+l)/R 
3/R; 00 

-(n+l)/(n-l) 
1; 3 

I -R/(n+2) 
0; R/4 

(9). 

(n+l)(n+2)/R 3 
12/R 2 ; “ 

(n+l)(n+2)/R(n-l) 
12/R; ~ 

-(n+2)/R 
4/R; ® 

1 



Table 1: Spectral formulas y for the operators connecting 

T, Ag, and , and their least upper 

^r ?>r 2 

bound (LUB) and greatest lower bound (GLB). 


hi Table 1 the quantities in the first column are expressed in terms of those of 
the first row, for example, 


Agn = 


R(n-l) /3 s t\ 
(n+l)(n+2) \^r 3 / n * 


The relations to the right of the diagonal show a damping of higher frequenc ies 
that indicates a smoother isosurface for the derived gravity quantity than for the 
original one. This principle is also expressed by the two bounds which have a 
finite range between zero and a certain maximum. The reverse rule is valid for 
the formulas to the left of the diagonal in the diagram. Here the higher frequencies 
are amplified and we speak of an unsmoothing procedure. Because the GLB are 
infinite the convergency properties of these quantities have to be carefully 
analyzed. A small amplitude in the spectrum for a frequency n- “ may increase 
to 00 by multiplying with an unsmoothing operator. Small errors in the original 
quantity may raise to a high amount in the derived data. 


Because the square of the spectral formulas of Table 1 
of the degree variances v_(T), v„(Ag), v n ^|T\, and v, 


Ej, and v n 0fTj 


shows the relations 
for the considered 
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gravity quantities, and for the observation errors too, the diagram is of additional 
use for the analysis of covariance formulas and error propagation. 


Assume, for example a gravity anomaly Ag is estimated with a certain 
amount of white noise; the higher frequencies of the noise will be damped in 
calculating the disturbing potential from the anomaly, the high frequencies of the 
noise will be amplified in a derivation. This is one of the well known reasons 


why the computation of 



from gravity anomalies is connected with more prob- 


lems than the geoid computation from Ag. 


For the discussion of satellite experiments the scheme of table one- valid 
for Surface quantities-has to be expanded in radial direction. The object will be 
the analysis of the height dependence of gravity quantities. We introduce three 
altitude levels , R near the earth's surface, r* and r 2 for low and high satellite 
altitude. The operator necessary for the generalization is the upward continua- 
tion operator, discussed later in detail. For the potential the spectral formula in 
spherical harmonics for this operator is (Heiskanen Moritz, p. 35; 33 (1-82 7 )) 


( 2 . 8 ) 


R(r 3 - R B ) 
A(P-P') 3 



y (B M (P / )R n .(P)+S B ,(P')S. B (P)) 

— J 

a— 0 


R . . . radius of an inner sphere (near the earth's surface) with surface a. 


r . . . radius of a concentric sphere with surface t pct, p'ea 
l . . . spatial distance between P and p ' . 

With the upward continuation operator of equ. (2. 8) derivation of the potential in 
altitude r from the potential given on a is possible. A more detailed discussion of 
the continuation of gravity quantities is given in Chapter three. For the spherical 
harmonic coefficients of T the upward continuation is carried out with equ. (2. 8) by-. 



where the indices R and r x of the coefficients indicate their relation to a and T. 
Similarly we have for the coefficients of the Legendre polynomial 

(2.10) Tl 1 T„ 

Thus, it is easy to generalize the diagram of Table 1 in radial direction, fix order 
to avoid confusion we omit the normal derivative — , which has about the same 
spectral features as Ag. dr 
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r 2 - altitude 
high 


r-i - altitude 
low 


R - altitude 
surface 



Table 2: Spectral formulas for the operators connecting T, Ag, and 

S 2 T and their [LUB;GLB] in three altitude levels R, r x , 
ars 

and r 2 . Only the operators in smoothing direction (damp- 
ing of high frequencies) are introduced; the direction is 
indicated by <3 ' . The spectral formulas in reverse 

direction where high frequencies are amplified, are reci- 
procal to these expressions; the LUB are inverse of the 
GLB in smoothing direction and the GLB are inverse of 
the LUB. 


To get an idea about the spectra of the quantities of Table 2, their fre- 
quency distribution from n=2 to n=70 for two different models for altitude 
h= 250km and altitude h= 0km are shown in Figure 1, 2, and 3, where h= r-R 
(R . . . mean radius of the earth). The two models are: 


Tscherning-Rapp (1974) 


var n (Ag)= s n+p ' 


425.28 (P-1) mgal 2 
(n-2)(n+24) 



Rbj e = 6369.8 km 
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and Kaula (1969) var n (C * , S *) = Y (C * 3 + S * 3 ) = 160 * a 10 

» = 0 11 


-12 


(*) dimensionless 


which becomes 

var n (N) = s 


n+l fkMl 3 160* 

L RG J n a 


10 -13 a R 3 

— — m* , s = 

r 


R= 6371km 


We assumed 


and 


Ag n = j r.m.s. (Ag) n j = | var n (Ag n )* | 
N n = |r. m. s. (N) n |= |var n (N)^|. 


By the formulas g inserted in Table 2 we computed for both models (T-R, and K) 
N n , Ag n , and f d T j , The result may give an impression about the spectra of 
Ur 3 J n 

these gravity quantities. The different results for the model Tscherning-Rapp 
(based on currently available gravity material) and the model of Kaula (based 

on the information available in (1966) especially for the component 

Sr 

indicate the problems that may be related in deriving an "a priori" analysis about 
instrumentation resolutions, error influence, and so on. 

Example: Derivation of T* from Ag„ x . 

A: Let us first look at the smoothing relation: 


Tn 1 = + ^Sn 1 ( bounds 0, r x ) 
then opposite to the smoothing direction 

with (a) and (b) follows 
t " = 


(a) 


(b) 


(c) 
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Figure 1: Amplitude spectrum of the Legendre polynomial coefficients N n 
of the undulation for the two models K and (T-R) and for altitude 0 and 250 km. 
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B: We first go in an unsmoothing direction: 




n+3 



then 


T n=-^i Ag n ;(0;R) 


The result from (a) and (b) will be the same 



(a)' 

<b)' 


(c/= (c) 


As Table 2 indicates the same problem may be solved by using several different 
ways, A and B are only the shortest of them. Without considering observation 
errors the result will always be the same. In addition, we see that at least one 
time an unsmoothing step has to be done in order to reduce the gravity quantity 
from altitude r x to R. As already mentioned this step has to be analyzed very 
carefully by considering the convergency properties of the involved data. 

The columns of Table 2 indicate the attenuation of the frequencies of sur- 
face gravity quantities by the factor n+1 , J n+S , or n+s , r > R, for up- 
ward continuation to altitude r. The higher frequencies are, therefore strongly 
damped with altitude. For that reason it is easier to get information about 
high frequencies from observations to a low orbit satellite than from a high alti- 
tude satellite. 

It becomes clear that under the assumption of comparable observation 
accuracy for the quantities involved in Table 2, the disturbing potential T R , 
which defines the geoid, may be optimally derived from R . For this 

second derivative no step in an unsmoothing direction has to be carried out by 
deriving T and disturbances in high frequencies are damped by the order of 
about 1 (exactly by jj , compare Table 1). On the other hand the high 
n (n+l)(n+2) 

frequency features in f ^ — T'| R will be very significant or what is the same it will 

L g r 2 J 

be strongly affected by local disturbances. An illustrative example for the dif- 
ficulties caused by these local disturbances is the trouble necessary for an 
accurate calculation of the topographic correction for torsion balance measure- 
ments. The derivation of a global data field with this instrumentation would be 
absolutely impossible. 

In order to get gravity information in short time for the entire world, 
satellite techniques are necessary. Their loss of information by attenuation 
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may only be compensated by a higher sensitivity in the instrumentation. 

The disturbing potential T and its gradient vt are not directly measured 
by satellite techniques but have to be derived from observations of the satellites 
position £ or velocity s . Therefore it is convenient to also introduce these two 
quantities into the spectral scheme of Table 2. The time dependent position vec- 
tor s and velocity vector _s are connected with the space dependent potential 
by the well known relation 

(2.11) 6s=VT 

with 6s . . . difference between the actual acceleration s and a ’'normal'’ 
acceleration s_ 0 computed from a low degree field. 

VT . . . gradient of the disturbing potential. 

The tracking elements to the satellite are direction, range, or range rate, i. e. 
components of the vectors £ and £ or linear combinations of single components 
of these vectors. 


Spectral analysis for the gravity quantities should therefore include the 
relations connecting s, s, ands. Kaula(1969) derived an interesting proportionality 
between the potential V, i.e. the space dependent integral of acceleration and 
the velocity s, the time dependent integral of acceleration. He gives? 


(2.12) v n (s) = 


GM ^ r a. 


r n(V*) 


where 

v n (s) . . . degree variance of the satellite velocity. 


Jb ... proportionality factor, G . . . gravitational constant, 
a* a 

M . . . earth's mass, a e ... equatorial radius, a . . . semi-major 
axis of the satellite orbit. 


v n (V*) . . . degree variance of the non-dimensionalized (*) surface 
gravity potential. 


The same relation holds true between the disturbing potential T* and the residual 
velocity 6s. The amplification of high frequencies in the spectrum of 6s in 
deriving 6s from 6s, will therefore be proportional to that for the derivation of 
VT from T. 

The spectral relations between s_, £ and is are in practice also dependent 
on the applied method of numerical or analytical integration and on the sample 
rate of the observations. It would be of great use to find closed expressions for 
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the amplification of the high frequencies in deriving £ from £ and £, for it^would 
enable us to estimate the error propagation from s or s to T£ , Ag n R , or 5 T R 

9r 2 n 

by inserting this propagation law into the spectral scheme of Table 2. 

Having equation (2.11) and (2.12) in mind the spectral formulas between 
T„ , Ag n and d T i n elevation r x or r B may be expanded by £, £ and s for 
5r 2 n 

which we don't have an explicit spectral expression for the connecting operators. 
This is shown in Table 3. 



r a - altitude 
high 




• • 
equivalent for the altitude r x (compare Table 2) 

Table 3: Interrelations between _s, £, and £ and their connection with 
the spectral formulas for the gravity quantities of Table 2 
in order, to get a more general impression of the smoothing 
and amplifying properties in the frequency domain. 


The position vector £ of a satellite is a very smooth function, its spectral content 
is concentrated at low frequencies. As we see from Table 3, £ will have about 
the same amplitude distribution in the spectrum as the disturbing potential, and 
the spectrum of the satellite acceleration £ is identical to that of the gradient 
VT by equ. (2. 11), and will therefore in the radial component be about the same 
as the frequency spectrum of Ag. It is therefore more difficult to derive the 
spherical harmonic coefficients for higher frequencies for Ag and also for T from 
position observations than from range rate measurements. From this point of 
view direct observation of £ would be optimal. 

Finally we try to reflect the character of satellite observations, which 
are carried out from one specific height level to another. Assuming the observa- 
tions have limited resolution but are not affected by noise, they may be seen as 
connections from the frequency level in elevation one to the frequency level in 


elevation two. In reality this means the measurements transfer the estimation 
accuracy from one altitude to the other. As we know, we need fewer coefficients 
to determine a satellite position with a certain precision for high than for low 
altitude. Thus, when a high satellite observes a low flying satellite by distance 
measurements with the postulated attributes, the position accuracy of the low satellite 
becomes the same as that of the high satellite. Because of this a certain amount of 
additional information about the higher frequency part is derived for the low alti- 
tude orbit. This principle is used in high-low satellite experiments and in 
altimetry. 

Summarizing the topic of this chapter we have derived a spatial spectral 
scheme (Table 1 and 2). Adding the smoothing properties of satellite position, 
velocity and acceleration (Table 3) allows us to analyze the essential attributes 
of terrestrial and satellite measurement techniques. The discussion is supported 
by also taking into account the characteristics of satellite measurements between 
different altitude levels. 

Because these considerations are only valid for a global representation 
of gravity quantities in spherical harmonics the conclusions cannot be applied 
without modification to methods of local improvement of the gravity knowledge. 

For this purpose a detailed analysis of the special technique in mind is necessary. 

3. Boundary value problem. 

In recent time, concepts for local improvement of our knowledge about 
the earth's gravity by satellite methods have become more important. If the 
benefit of the higher accuracy of new techniques is utilized any global solution 
is coupled with a rapidly increasing number of unknowns in a least square adjust- 
ment. For a representation of the gravity field by spherical harmonics where 
the coefficients are the unknown parameters, the number of unknowns for a 
representation up to degree and order 36 (= 5° equal area solution) is already 
1363. For this reason a local adjustment where only mean values for a limited 
number of surface blocks are derived would be preferable. This kind of 
solution seems to be in contrast to some theoretical considerations which have 
their origin in the downward continuation problem for gravity quantities. 

For a better understanding, a short review of the boundary value problem 
will be given: 

Problem statement: Given any function w oh a surface S, derive a function f 
defined outside S and harmonic which approaches on S the given boundary function 

w. 

This first boundary value problem or Dirichlet's problem has an easy 
solution for a boundary sphere a. (R ... radius of the sphere, 2 . . . space 
outside 0) in form of the Poisson integral, 
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(3. 1) 


f(P) = 


R(r 2 -R 2 

4rr 



a 


UP') 

l( P, P' ) 3 


da 


(All terms except f are already explained for equ. (2.8), where the Poisson 
operator was introduced into the spectral diagram). 


The Poisson integral solves the Dirichlet problem for the sphere a; that is, 
given any continuous function f(P 7 ), P'ea, the Poisson integral of f defines a 
function exterior to the sphere, and harmonic in 2 which approaches the given 
boundary function. 


Given the boundary function f(P 7 ) on a the integral formula (3. 1) allows 
the determination of f(P) for all P££ . In this formulation equation (3. 1) defines 
the upward continuation problem. 

But far more complicated is the inverse problem where our aim is the 
derivation of the boundary function f(P') from f(P). Equ. (3. 1) becomes an 
integral equation of the first land and the problem is named improperly posed. 
For a general but given boundary surface the problem was analytically handled 
for example by Week (1972) and Krarup (1969). Discrete and statistical solutions 
of this so called downward continuation problem are proposed by B jerhammar 
(1966), (1969), Moritz (1970) or Schwarz (1971). The solution of the integral 
equation is essentially facilitated by choosing a spherical boundary surface a as 
we already assumed for the Poisson integral of equ. (3.1). Equation (3. 1) may 
be more abstractly written as: 

(3.2) f(P) =j^f a(|P-P'|> f(P # ) da 

a 

with 

a(|p-P r |) ... Poisson operator, where | P-P* | indicates that a(x) 

is a function of only the distance |P-P / j. 

Expression (3.2) formulates a convolution integral, it has the advantageous 
property that it becomes a simple product after transformation into the fre- 
quency domain. Functions defined on the real line are transformed into the 
frequency domain by Fourier expansion. The analogic procedure for a function 
defined on a sphere is the expansion into surface spherical harmonics. 

Our convolution integral will therefore become by inserting the spherical 
harmonic expressions for f(P), f(P') and a( | P-P' | ) 

(3.3) f M = a nm C 
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Recalling the spherical harmonic expression for the Poisson operator (2. 8) 


( 2 . 8 ) 



<R nB (P')R nB (P)+S nB (P')S nB (P)) 


and expanding the potential V into spherical harmonics 

00 n 

(3.4) V(P')= £ £ (C(V)nm Rn B (P') + S(V) n R B S nB (p')) 

n— 0 m=0 

we obtain by inserting equ. (2. 8) and (3. 4) into equ. (3. 1) with V(P) for f(P) and 
V(P') for f(P') 



R 


n + l 


C(V)n # 



\s(V) R J 



which is in accordance with equ. (3.3) and the expression for the disturbing 
potential T, equ. (2.9). Also the inverse formula - the downward continuation - 
holds 



In the spectral representation an easy solution of the downward continu- 
ation problem in the form of equ. (3. 6) seems to exist. But in order to be able 
to expand a function f (P ) into spherical harmonics it has to be given, independent 
of the frequency (n, m) in which we are interested, as a continuous function cover- 
ing a sphere, expressed by 


(3.7) 


c <f)n„=^ /f(P')R nB (P , )do 

i a _ 

S(f)n „=^ f f(P')S nB (P')dcr 


a 


Because the coefficients derived by equ. (3.7) are independent, one from another, 
it is possible first to estimate the coefficients up to a low degree and order and 
then to complete them step by step. Thus, when we have given a smoothed 
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form of f(P / ) that contains only frequencies up to degree un we may derive the 
coefficients c(f) and s(f) up to degree n« by (3.7).* For such a procedure a solu- 
tion of upward and downward continuation of gravity information in form of 
spherical harmonic coefficients is possible. These theoretical considerations 
form the background for actual global satellite solutions. 


But there is a gap between theory and practice: The data from which the 

coefficients are derived does not cover a global sphere with a certain radius but 
vary in altitude because of the elliptic shape of the orbit; will not cover all sub- 
regions with the same density; and is given in discrete but not in continuous form. 
Thus the integral formulas (3.7) have to be approximated by a system of dis- 
crete linear equations applicable for least square adjustment. The coefficients 
for the spherical harmonic expansion will therefore lose their independence 
and will be affected by systematic errors. In addition,not incorporated higher 
frequencies in the data will falsify low frequencies by alaising which also raises 
the correlation between different coefficients. 


More difficult is the problem of local improvement of surface gravity 
information where we are directly faced with integral equation (3. 1). Frequency 
representation can only be used in a modified form. 


4. Satellite to satellite tracking and satellite gradiometry. 


The large number of unknowns connected with a refinement of our know- 
ledge of the gravity field by global analysis and the varying density of satellite 
observations depending on the orbit elements, the technical equipment, and the 
distribution of ground tracking stations leads our attention to local methods. 

This idea is supported by the fact that terrestrial measurements are given with 
high accuracy and density in some areas with no data in other areas. Our 
main interest lies in gravity anomalies from which by combination with ter- 
restrial anomalies a geoid computation can be carried out. The satellite techniques 
in mind for this purpose are satellite to satellite tracking (SST) and satellite 
gradiometry . 

We do not include satellite altimetry into these considerations, though 
this technique is very promising. But the altimeter observation connects 
directly the satellite orbit in altitude r with the topography of the earth, i. e. 
with respect to the discussion of chapter two, downward continuation from alti- 
tude r to the surface is carried out by the observation itself whereas we want to 
analyze the problem connected with the analytical downward continuation from 
satellite altitude. 


Directional and range measurements from ground tracking stations are 
not accurate enough to provide boundary gravity anomalies for smaller block 
size with sufficient precision. The high frequencies in the spectrum of the 
disturbing potential T in satellite altitude rare damped by the factor £Rj n+1 
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as shown in Table 2. By deriving the satellite position s_ from the velocity £ 
which is proportional to T, equ. (2.12), the spectrum is again damped. The 
instrumentation for direction and range measurements which are components of 
the satellite position, are therefore usually only sensitive enough to resolve the 
low frequencies of the disturbing potential. 

Studies by Forward (1971) and by Glaser and Sherry (1972) assume the 
same frequency properties for SST and the disturbing potential and also for 

gradiometry and d T . Therefore we may expect a low frequency improvement 
?ir 8 

from SST and better knowledge of higher frequencies by gradiometry. These 
conclusions are based on the essential simplification that SST has the same fre- 
quency features as s itself and the only important component for the disturbing 

potential tensor is the radial component f. Z . 

dr 8 

When we try to avoid any misleading simplification we are faced with the 
complex formula systems connecting the observation data with the final gravity 
anomaly. In spite of that we will try to obtain an insight into the problem by 
dividing it into several steps and analyzing their individual features. 


4. 1 Satellite to satellite tracking 

The problem with which we are mainly concerned here is to derive gravity 
anomalies Ag(p') from a special type of satellite observation O(P), P'coand Pct. 
The most common way of establishing a linearized system of observation equations 
is, Kaula (1966), 




C W 

n + nx s 


dS 


O . .. n-dimensional vector of observation 
dO ... " " " residuals 

£ . . . " " " the calculated observation elements 

W . . . nxm sensitivity matrix (coefficient matrix) 
d.$ . . . m-dimensional vector of the unknown parameters 

For range rate measurements R equ. (4. 1) will become 
(4. 2) R 0 _ c + dR « W Ag , with Ro-c = R 0 - R c and 


Ag . . . vector of unknown gravity anomalies. 

The sensitivity matrix Wconsists of the coefficients 


5R t 

9Ag 3 


and shows the 
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interrelation between the range rate observation at a certain moment and a 
distinct gravity anomaly. As we know, equ. (4. 2) can be divided into different 
parts by expanding 

(4.3a) R 0 .c+ dR = Ado; 


where A is the sensitivity matrix with respect to the intermedium parameters 
da and 

(4. 3b) da= B Ag 

where the coefficient matrix contains the relations between da and Ag. By com- 
paring (4. 2) and (4. 3 a,b) we see 


(4.3c) W= AB. 

Applied to our problem we may first analyze the evaluation of the acceleration 
s(P) or what is the same VT(P),( compare equ. (2. 11)), from range rate data 
R(P). Then we try to clarify the derivation of Ag(P 7 ) from VT(P) . For a term 
w lk with indices (i,k) of the sensitivity matrix W for range and range rate ob- 
servation in SST we get, compare Martin (1972) and Hajela (1974) 


, . m4“> 8tV 8Ri HU> 8T„ r ' J 8R< 3u) 9pV 

(4.4) w lk = + — — and 

ST,' 1 , Mg k ^T x r \, 3Ag k BT x r s a 3Ag k 


• fo) r. J . feu) r, J . feu) r m j 

? . 3Ri 3T X X 3Rj 3Tx 3R t ■ 3T X * 

(4.5) w lk = “ + — ~ — 

*T, X , 3Ag k 3T x 1 < 3Ag k 3T x ? j *Ag k 


with the index conventions 


l . 

k 

3 


J 


range or range rate observation number i with respect to 
gravity anomaly number k 

sub and superindex j means summation for j = 1, 2, 3 of 


the vector components T, 


T x- , and T Xa of VT 


a 3 - * ... 

r x , r a . . . as already introduced, height level r x (low) and r a (high). 

Id . . . observation from relay satellite down to the close satellite 
2u . . . observation from the close satellite upward to the relay satellite, 
compare again Hajela (1974). 


Each term 


3R 

3T X 


is formed by 


19 


3R = _^R_ SX A _3R_ Sxf 

( 4 * 6 ) 3T X 9X A dT x aXj aT x 

and each term -||r by 

,4 71 22 .is. 

' ' ' ST, 'X £ 8T, 


where 

£ 


is summation for £=1,2, 3. 


According to equations (4. 3 a-c) we partition the sensitivity matrix W 
into a product of three matrices A,B, and CX Matrix A contains the elements 
dll and _®_for range rate measurements, equ. (4.6), and the elements dR 


ax, 


t 


ax„ 


ax« 


are 


for^range observations, equ. (4.7). Analogously the coefficients of B 

5X and dX tV»A ranore rflfp nnH ax x fnt' flip nrTO m TVin fow 

aT x aT x 

STv 1 

and . We divide the sensitivity analysis for Ag estimation from range 


for the range rate and - dX x for the range method. The terms of matrix C are 

aT x “ 




dAg k aAg k ‘ 

or range rate observations into the analysis of the matrices A, B, and C.. 

dR 3R BR 

Matrix A: The coefficients , and ^r- are the projections of the range 

£ £ £ 


or range rate vector onto the components of the position and velocity vector of 
the satellite, expressed in the inertial coordinate system (Xi, Xg, X3). Their 
magnitude and therefore their influence on the final result is dependent on the 
geometrical configuration of the two satellites. The sensitivity coefficients in 
A may be transformed either into an earth fixed (<p, X, r)- system into radial, 
longitudinal, or latitudinal direction or into a radial, along track and cross track 
component. In the latter system Koch and Argentiero (1974) carried out sensi- 
tivity considerations for the GRAVSAT/GEQPAUSE mission. 

tin 

In a high low SST experiment the projection onto the radial coordinate axis 
of the range or range rate vector between two satellites will be predominant. 
Therefore the radial variations in range and range rate will be essential for the 
gravity anomaly determination. 

The geometrical configuration in a low -low SST experiment will change 
continuously, no predominant coefficient occurs. The three components with 
respect to the axis of a certain coordinate system will contribute equally to the 
final result regardless of their importance with respect to Ag. 

Equations (4. 6) and (4. 7) show the presence of derivatives with respect 
to the satellite position X and the satellite velocity X for range rate experiments 
whereas for range measurements, as expected, only derivatives with respect to 
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X are involved. This is the essential difference in the solution equations for 
both types of observations. As discussed in detail in Chapter 2, the satellite 
position is only in a very smoothed form dependent on the gravity field. For 
example, Chovitz (1973) derived that the effect of high harmonics in the poten- 
tial (n£ 4.0) on a satellite position in altitude 900km is absolutely negligible. 
Sufficient sensitivity for the gravity anomaly determination may be only expected 
from satellite velocity dependent terms. Therefore in the following cur consid- 
erations are limited to range rate experiments. 

Matrix B: The coefficients — — and in B connect the satellite 

! &T XJ ST XJ — 

position and velocity vector with the components of the gradient vt of the 
disturbing potential. The gradient VT is equal to the acceleration vector X. 

The magnitude of these terms is dependent on the method of integration from 
Xt to X t and X t at every instant t. For high precision in a numerical integration 
many terms of type X t+1 , i= -k(0)A, have to be introduced into a linear combin- 
ation. In order to avoid high error correlation between these neighboring data 
and a considerable loss of resolution caused by the method of numerical inte- 
gration a high sample rate is desired which has its limit in the technical equip- 
ment. 


9T^ 


5T? 


The transformation of the 


Matrix C : The elements of C are u ^ and 

3A Sk Ag k 

components T x , j= 1,2,3 of VT in the inertial system to the components 

3 bt w 


ajT » , and T^= ^ in the earth fixed (r, <p, X) system is expressed, 


Rapp (1974, p. 5), by 


f T * ' 


■*3r acp *X ’ 


T r ' 

1 


^Xx BX x ax x 



. 

' 


ar ax 



t ’ 2 


?)Xg ^Xg 9X7 


T 

a <£> 

i 

rp 


ar a<p ax 



1 v 

X ^ 

L J 


L9 x 3 ax 3 ax 3 J 


T x 


The gradient in spherical coordinates _VT, Ph has the form 




T r 

— T 
r <P 


_1 

ycos <p 
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and is related to the gravity anomalies by the expression,) Heiskanen-Moritz, 
p. 234 | ; 




— ““Ii 


(4.8) 


- f [ Ag — cos a da 
r J J Kip 

a 


with St(r,j/)) . . . Pizzetti's extension of Stokes' function, (Heiskanen-Moritz, 
p. 92, Hotine, 1969, p. 310) 


It is easy to find from (4. 8) the derivative of T x with respect to a cer- 
tain mean anomaly Ag k by approximating the integral by a finite sum,(Hajela 
1974). The influence of each of the three terms of equ. (4. 8) in their linear 
combination expressed by C is decided by the magnitude of the derivatives of 
Pizzetti's function with respect to ip and r and the area of the finite surface 
elements Aa for the block mean Ag k . The term dSt(r , ip) i n the T<p and T^ component 

is rotation invariant and projected by cos a and sin a onto a desired azimuth. 
Therefore a comparison of - — dStjr, ip) and dSt(r,4)) sufficient. 

By considering the spectral properties of these two functions we get an 
impression of their behavior in the frequency domain. Both functions are con- 
nected with Ag by the convolution integrals in (4. 8) and therefore act like 
filters which decide about the frequency content of the derived gravity anomaly. 

Pizzetti's function St(r, tp) is expressed in an expansion into Legendre 
functions, by (Hotine, 1969 , p. 311) 


00 

(4. 9) St(k , ip) = Y — k n+1 (2n+l) P B (cos ip) 
Li n-1 


n=3 


R 

which becomes with k= ~ and k<l 

' ' ' . '.00 

(4.10) St(r, ip)=B y - 4 - 
Aj n-1 


n= Z 


u+1 


(2n+l) P n (coSi/j) 
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The radial derivative is 


(4.11) 5^3 
dr 


00 




n + 2 

(2n+l) P n (cos ip) 


and the derivative with respect to if) 


(4.12) 


1 mjrjJ)) 

r bif} 


00 



n= Z 


1 R 
n-1 r 

M J 


n + a 

(2n+l) 


dPJcos ([)) 

dlf) 


In the radial component of VT EpU the spectral content is connected with the spec- 
trum of a certain gravity anomaly Ag by the factor Bii. f il n+s . Thus, be- 

fro in +8 k n-1 


side the attenuation effect 
because lim = 1. 

n -> 00 II — 1 


(?) 


the two spectra are in a one to one correspondence, 


For both other components of VT , which have only poor high frequency infor- 
mation, the derived amplitudes in the Ag spectrum are amplified by (n-1), because 
of the factor in equ. (4. 12). It will be difficult to obtain from these two terms 

higher frequency information in Ag which is not seriously affected by errors. 

Since the high frequencies at high altitudes, where the relay satellite is 
located are very strongly damped , the ST x r 3 coefficients will hardly contri- 

SAg k 

buteto the derivation of gravity anomalies of small block size. The influence of 
each of the components of equation (4.9) is dependent on its magnitude in the 
linear combination which is expressed by the matrix product 


WAg = ABC Ag 

and is by A mainly a function of the geometrical configuration of the two satellites. 
4. 2 satellite gradiometry 

In satellite gradiometry the equations are more transparent. The gravity 
quantity is measured directly by a satellite mounted instrumentation. Whereas in 
SST the gravity gradient was derived from range rate perturbations on a predeter- 
mined satellite orbit. 

Two different gradiometer types are in discussion, the hardmounted and 
the rotating gradiometer, (Forward, 1971, Reed, 1973 ) . For the amplitude of 
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the output of both types an extended sensitivity analysis was done by Glaser 
and Sherry (1972) and in more detail by Reed (1973). 

We are interested in the resolution of the measurement elements with 
respect to the spectral content of gravity anomalies. The components of the 
gravity tensor are measured in a moving frame. By coordinate transformation 
they are transformed to the earth fixed (<p, X, r)- system. Thus, in a least 
squares solution the sensitivity matrix W will again consist of several terms 
with more or less desired properties. 


The interrelations of the tensor components of T in the local rectangular co 
ordinate system (7), £, C) — K ... radial direction, r\ . . . tangential latitudinal 
direction, %... tangential longitudinal direction to the components in the earth 
fixed (<3, X, r)- system are derived in Reed (1973), 


T 35- 

l rm 


COS (0 


1 _ v tan<p 1 _ 1 

< 7 T xx> + ~r <; T <p ,+ ; 


t « = <7W + T t ' 


T ? £=T rr 


(4.13) 


sinto 


T £f? T r?£ cos<p ^ ~P V + rcos a <p ^r T X^ 


T t?C cos <p ^r T Xr ^ rcos<p ^r T X^ 


T„. = T>v> = ( — T 

££ r <pr 


)- V 

r r 


The connection of the second derivatives of the disturbing potential in the (<p,X, 
r) -system to the desired gravity anomalies are again expressed by derivations 
of the extended Stokes' formula with respect to r, <pand X, for example 


a 

Therefore the discussion of the C -matrix for the range rate experiments can be 
repeated analogously for gradiometry. To get an insight into the frequency 
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behavior of the kernel functions — ^ (representing the kernel 

functions for T r(a and T r \ ) and - (representing the kernel functions for 

v dip 3 

T , T. . , and T. ) we derive their spectral formulas: 

<po XX Xtp 


5 3 St(r, 0) = 1 r 1 (n+l)(n+2) 

•s__3 / > 




dr 3 

r -■ (n-1) 

n=s 

(4.14) 

1 

5 2 St(r, 0) 

= _ 1 y n+1 

(a-c) 

r 

3r 50 

r L n-1 




n= 2 


1 

a 2 St(r ,0) 

00 

= 1 
r L n-1 

n = 2 


r 3 

30 s 


n + 2 


(2n+l)P n (cos 0) 


v “* = 

r 


(2n+l) 


3P n (cos 0) 

dip 


R 


n + 3 


3 2 P„ (cos 0) 

( 2n+1 )~ s? — 


As already discussed in Chapter 2, the second radial derivative has very good 
high frequency information. Errors in the radial tensor component T rr are 
damped in deriving Ag by the factor ( n ^* n+2 j> ec t u - ( 4 - 14a), and high resolution 

can be expected. For the purpose of Ag - computation in small block sizes (which 
have high frequency content), the T rr component seems to be well suited. 

All derivatives of T in the tensor with respect to r and an additional com- 
ponent (<p or X) are in a one to one correspondence with the gravity anomalies, 
compare (4. 14b). 


For second derivatives not containing a radial component, the kernel (in 
equ. (1. 14c)) function will amplify the frequencies of the tensor components in 
T(p<p, and by (n-1). These tensor components have only good low fre- 
quency information. Therefore we may expect an error influence on the higher 
frequencies in the mean block anomalies for these components. 


The other members in (4.13), T r , ~ T^ and i.T^ are all damped by ^ . 
Their influence on the resulting gravity anomalies will be small. 

The signal amplitude for the rotating gradiometer is, Reed (1973), 

| 

amp= ((V 3 3-V 11 ) s +4V 1 2 3 ) , 

subindices 1>Sf3 ... moving frame, 1-axis directed outward. 
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The component V u would have the advantageous properties of T rr but cannot 
be separately evaluated. By the combination with V 33 and V i3 the result is de- 
pendent on the worse properties of these terms. The deduced gravity anomalies 
cannot be resolved with the same accuracy than from T rr alone. 

Only for the hard-mounted system it seems to be possible to separate the 
radial derivative T rr if an accurate space orientation is available. In a closed 
adjustment for five independent tensor components the special features of T rP 
would be disturbed. But as far as we can see only the rotating system has a 
chance to be realized in the near future with sufficient accuracy. 

These considerations show that the advantageous properties of T ( p ) an( j 

or 2 

for a high resolution in gravity anomalies, indicated by the discussion of 
or 

Chapter 2, are with the up to now available techniques not realizable. Sensitivity 
studies based only on the radial components are therefore to a certain amount 
misleading. 

The aim of the investigations in this chapter was to demonstrate the 
advantage of an "a priori" - sensitivity analysis of an experiment. The coef- 
ficients of the sensitivity matrix W can be evaluated by the mathematical inter- 
relations between measurement and unknown quantity. No simulations or practical 
data are assumed. For SST it seems to be possible to deduce objective criteria 
about the optimal geometrical configuration of the relay and the tracked satellite, 
a possible need for more than one relay satellite to get a better geometrical con- 
figuration and the requirements in the numerical integration. An accurate compar- 
ison of the sensitivity matrices in SST and gradiometry should give an impression 
of the individual advantages of both techniques. Detailed conclusions may be drawn, 
as indicated by a similar "a posteriori" sensitivity study for traditional satellite 
techniques by Gaposchkin (1970). 

5. Downward Continuation from Satellite Altitudes. 

As already mentioned in Chapter 3 we are mainly interested in deriving 
surface gravity anomalies from satellite measurements. In order to limit the 
number of unknowns a local improvement seems to be more promising. The 
derivation shall be done from the gradient of the disturbing potential VT by 
satellite to satellite tracking techniques or from the components of the anomalous 
potential tensor in satellite gradiometry. 

Up to now this topic received little attention as a downward continuation 
problem. Theoretical and numerical studies treated the continuation from 
aerial altitudes or the reduction from the earth's surface to the geoid, for 
example Dean (1958), Madkour (1966), Koch, K-R. (1968), Moritz (1970), 

Schwarz, K-P. (1971), Bjerhammar (1966) and (1973). 
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An impression of the downward continuation effect on the gravity anom- 
aly variance for ranges of n which are essential for 5 u x5°, 2.5°x2. 5° , and 
l°xl° anomalies is given in Table 4. 


The formula used for Table 4 and Figure 4 is: 
var(Ag) Si ,B, = ^ 


2i>+4 


A 

n+B 


(A) 


which is a simplified version of the expression 


"a 


var(Ag) ni , n = ) s 


n+2 


A(n-1) = 

(n-2)(n+B) ’ r 3 


(B) 


derived by Tscheming and Rapp (1974). It may be evaluated either computa- 
tionally in a summation routine or analytically. Then we get 


var(Ag) n no = y s n+3 & ± As 2 f dn 
1 ^ n+B J 


n+B 


= As a [e" B ^ 18 Ei (0ns (n+B))] 3 , 


because / ■ = £ e~"~' ~ Ei ( - (a+bx)), (Meyer Zur Capellen W. , Integral tafeln, 

■' a+bx b b Berlin, 1950). 


k kx 


1 „-*»/ b ;rr, k 


The exponential integral Ei(x) is listed in standard tables like 
(Akademija Nauk SSR, Tafelnder Exponentialintegrale, Moskau, 1954). 

Two facts are of particular interest: 


1. The changing ratio of the high to the low frequency part for different 
altitudes, for example 


var ( A gk.8 Q - ^ whereas var ( A S> 3 -s° = M 

var(Ag) 73fl80 h _ okm 1 var(Ag) 73flB0 h=2 00km 1 


2. The immense attenuation effect for higher frequencies as can be 
seen in Fig. 4. For example the amount var(Ag) 73iieo i.e. the additional infor- 
mation content distinguishing a l°x 1° anomaly from a 2. 5°x 2. 5° anomaly 
(assuming the simplified rule that the highest frequency degree n in a 0 x 0 
anomaly is n= ) is amplified from 0.6 mgal 2 at 200km elevation by 510% 
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Table 4 


Sum of Degree Variances for Gravity Anomalies var(Ag) v Between 
Degree n x and tu, for different altitudes. 


\ [mgal 3 ] 

[km] \ n = 

altitude's^ 3 2Q 36 72 180 

0 

235.7 

129.2 

194. 5 

305.2 


100 

163.9 

51.1 

36.8 

10.3 


200 

117.8 

20.9 

7.8 

0.6 


300 

87.1 

8.8 

1.8 

o.o 4 


400 

66.1 

3.8 

0.5 

0. 0q 00 4 




Figure 4: Attenuation of the variance in Ag between 
degree n x and degree n, from altitude Okm to 400 km, 
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to 305. 2 mgal 2 at the earth's surface. An observation error of only ±0. 5 mgal 
in this frequency range (at 200 km) will rise to 128 mgal 2 at the surface and 
cause an uncertainty of about 15% in the l°x l u anomaly with respect to surface 
variance of 865 mgal 2 . These few numbers illustrate what extreme accuracy 
should be reached by the observation equipment for a l°x 1° anomaly resolution. 
It also explains why it is necessary to discuss the topic as a downward contin- 
uation problem. 

The center of the following considerations will be 


(5.1) T(P) 


4tt / St <r ’ 


Ip) Ag(P')dor , 


with St(r ,lp) . . . pizzetti's extension of Stokes' function 
p'ca ... as in Chapter 3. 

Equation (5. 1) connects the anomalous potential T in a space point P€S with the 
boundary function Ag(P') on the sphere <7. In what way the continuation oper- 
ator of equ. (2.9) is involved becomes obvious by expanding (5.1) into spherical 
harmonics. From (4. 10) we get for St(r,$), 


CO 

(5.2) St(r, 0)= V -±- 
Lj n-i 


n“ s- 


n + l n 


) (®ne (P)Rn« (P*) - ** Sn* (P)S nB (P # )) 

n~ 0 


and we end up in analogy to equ. (2.6) with 

^ c (T)nn 
.s(T) na 


(5. 3) 


>= R 


n-1 


n+l j c' (Ag)„ B ^ 

Vg)n B J 


The same result can immediately be derived from Table 2 by connecting Ag nn 
with T„^ . The inverse formula is 


(5.4) 



r2i n_1 > 

V 

r 

r 

R 



The problem is to solve the integral equation of the first kind (5. 1) for a finite 
number of gravity anomalies given discrete measurements of the gradient or 
the tensor components of T. Integral equation (5. 1) formulates an "improperly 
posed problem", as the inverse operator of St(r,j|)) is not bounded as may be 
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seen from (5. 4). For an analytic study methods of an approximate solution are 
indicated by Week (1972) and Bjerhammar (1973) or deduced in more detail for 
the planar approximation by Schwarz (1971). On the other hand there is also 
proposed the method of collocation in the statistical or analytical sense (re- 
producing kernel Hilbert spaces) by Moritz (1970) and Krarup (1969). 

As we will see there is a lot of complications connected with the deter- 
ministic approach to find a solution of the downward continuation problem by 
means of equ. (5. 1) whereas the solution by prediction methods seems to be 
very easy and extremely elegant. Logically one may ask why the prediction 
shall be so much easier. Where have the difficulties for the deterministic 
solution gone in the prediction approach? 

Krarup (1969) defined a Hilbert space H of all given potentials <p regular 
in 2 with square integrable boundary values on a. (This last restriction is 
essential and enables a solution of the downward-continuation for <p to be 
found). In addition, with the assumption of rotational invariance an expression 
for the covariance function of the disturbing potential is: 


P n (COSj/)pQ ) 

where P and Q may be both €2 or ea, 

v n (T) . . . degree variances of T (units r potential) 

The kernels in the integral formulas of physical geodesy are usually known 
from the solution of the boundary value problems, hi a first approximation 
these boundary value problems are solved for a spherical boundary surface. 
The spherical reference! surface generates rotational invariant integral ker- 
nels with a reproducing property in H. The reproducing property of the 
operators allows us to der ive auto-and cross-covariances for all other quan- 
tities related to the earth’s potential. The Poisson kernel, eq. (2.8), to be 
applicable for both P and Qc2 is generalized to 

n + 1 

(2n+l)P„ (cosi/v Q ) 


p f R 2 

(5.6) Kpdi(p, q> - r ^ 


(5.5) Ctt (PjQ)= ) W) 


R 

r P r Q 



n- o 


R c 


n+ 1 


I r p iqJ 


X (Rn n (P)Rn« (P / )+S nu (P)S no (P')) 
0 


A linear prediction approach for the derivation of gravity anomalies 
Ag from disturbing potential values T may be formulated neglecting measure- 
ment errors by 
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(5.7) Ag P / - B T P 

Agp ' . . . vector of unknown quantities 
B . . . unknown system matrix 
Tp ... given data set 

Given Tcortinuous on a sphere T with radius r in satellite altitude and 
concentric with ct equation (5. 7) becomes 


dT 
T 

First we try to find an expression for the system function B(P'-P). In a sta- 
tistical approach we take the expectation E(.) on both sides of (5.8) with 
respect to the given function T(Q), Qct. 


(5.8) Ag(P / ) = J B(P'-P) T(P) 


(5.9) E(Ag(P / )T(Q))=y' B(I^-P)E(T(P)T(Q)) dT 


Assuming rotational invariance for the "stochastic’' process T(Q) we insert 
Ctt from equ. (5. 5) for E(T(P)T(Q)). For E(Ag(P')T(Q)) we introduce 


(5.9a) C AgT (P',Q) 


00 

> 


- (n-l)v„(T) 

rC 


R 

r Q 


n + 1 


P n (COS&' Q ) 


The unknown system function related to a sphere can be developed into 
Legendre polynomials 


(5.9b) B(P / -P) =2 J x n P n (cos0 P 'p) 

V n= 0 

Inserting (5. 5), (5. 9a), and (5.9b) into equation (5.9) gives using the ortho- 
gonality relationships 
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The x, terms are exactly the coefficients of the spectral representation of 
the "downward continuation" operator we derived in equ. (5.4), We see that 
for the global case the solution of the linear prediction — assuming no data 
errors — has the same result as the analytic approach. 

The basic equation is the "Wiener-Hopf" equation (5.9) which may be 
rewritten for this example as 


( 5. 11) C Ag f (P' , Q) = -P)Cr T (P, Q)dr . 

T 

By (5. 11) the previously mentioned reproducing property of the integral kernels 
in the spherical approximation of the basic equations in physical geodesy be- 
comes obvious. Equ. (5.11) is the expression for the so called "covariance 
law of propagation" which allows us to derive the auto-and cross covariances 
of all harmonic quantities related to the earth's potential from C T r of equ, (5.5). 


Because of the orthogonality relationships, equations (5.8) to (5.11) hold 
true also when they involve functions containing only frequencies up to a finite 
degree N, as it is applied for global satellite solutions, compare with equ. (3. 7). 
Therefore, the prediction method will not beconnectedwithanessential advantage for 
global appl ications compared w ith the determ inis tic method . But already for combina- 
tion solutions with terrestrial data prediction maybe more convenient. For a finite 
amount of data not covering the total sphere the integral formulas bee ome matrix equations 
which approximate the continuous expressions only to a certain degree. 


The "Wiener-Hopf" integral equation degenerates to the matrix equation 


(5-12) C = 2 C , p , 


and 


P 'Q 


-1 


(5.13) B C_Ag P ' t q £ t p t ( 


When we insert the finite expressions of Ctt and C A g T into equ. (5. 12) 
the orthogonality relations are lost and B will not be the finite, discretisized 
form of the "downward continuation operator". Logically the deterministic 
approach, derived from the discretisation of equ. (5.1) and the prediction approach, 
deduced by inserting equ. (5. 13) into equ. (5. 7) will lead to different results. In a 
smoothing process such as upward continuation or calculation of undulations from 
gravity anomalies the continuous formulas are integrals , not integral equations. 

Here the discretisation cannot be very misleading. The discretisation error will 
be about the same magnitude as the prediction error. Generally, prediction will 
result in smoother quantities — in this case a desired feature — and can be easily 
carried out for non-homogenuous data distribution. 

The solution of the downward continuation problem is an unsmoothing pro- 
cedure in the manner of the computation of gravity anomalies from undulations. 

For such problems the approximate solution of an integral equation is desired. 


22 


Equ. (5. 1) becomes: 


(5.14) T p + n = A Ag P / 

n, . . . observation noise, with E(n) = 0 and E(n n' ) = C nn and the 
problem is to find a proper inverse of A . A conceptual analogy between the 
question of finding an optimal inverse of_A by linear estimation in terms of 
generalized inverses and the collocation (generalized prediction and adjust- 
ment) concept was noted by Moritz (1973). For problems of this type the pre- 
diction method seems to be more elegant because it apparently avoids the 
trouble connected with finding the proper inverse of A. 


5. 1 Least Squares Prediction Solution 

The prediction formula for the discrete downward continuation will be 
(with equ. (5. 7) and equ. (5.13)): 

(5.15) I|,' =£a & ',, C V q T, 

By taking into account observation noise n as in equ. (5. 14) with E(Tpn')=0 
equ. (5. 15) has to be modified to 

(5. 16) Ag P 7 = C ^ C Xp and Xp = T p + n p 

£ = £ t p t q + £»p» q 

For all practical applications T P will consist of a very limited number of ob- 
servations. Now CT p Tq expresses the real physical relation between the ob- 
servations and C_^g T the real physical relation between the observations and 
the desired unknowns' 5 . An autocovariance matrix "C Tp T p estimated from the 
observed values T^ P may not have a high confidence because of the limited 
number of samples and C^g p/T cannot be calculated on the basis of the given 
datas for the Ag P are unknown.^ Therefore the sample covariaices are 
supplied under the assumption of isotropy and stationarity by the global covar- 
iances. In addition these general covariances make the results of the limited 
number of T P consistent with later observations in other areas and a global 
solution. 

Statistically speaking a continuous correlated process on the sphere is 
approximated by a discrete n-order Markov model. For this model the optimal 
solution would be the conditional expectation 
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4p' = E(Agp' | Tp+ n P ) 
with the conditions 

A 

E(Ag|/) : = E(B(T P +n P ))=E (Agp / ) [unbiased] 


and 


EjAgp' - Agp'| * ^ E|Agp'_ B(T P + np)|^ 

with basis C . . . nonnegative definite and symmetrical. 

In our case the first condition is not fulfilled because Ag is "absolute" 
defined on the basis of T(P). Therefore it is not desirable to solve the prob- 
lem on the basis of sample covariances which would lead to a biased result. 

The drawback of introducing the global covariances into equ. (5. 15) is that 
their system of eigenfunctions and eigenvalues is not the system of eigenvectors 
and eigenvalues of the discrete finite solution. The results will be very smooth 
and optimally consistent. Smoothness is not at all desirable in downward con- 
tinuation where we try to find the unsmoothed quantity from the damped. But 
the control of high frequencies by prediction will bound the errors and there- 
fore make a solution of the "improperly posed problem" feasible. The follow- 
ing conceptual disadvantages are essential for the prediction application. 

— The finite approach predicts Ag P only with the least possible pre- 
diction error €p= ( Ag , - Ag,) dependent on the introduced covariances 

(5. 17) var(Cp ) - C A g p/ A g p ' - C^/ Tp Ct p t„ C ^ 

The variance of the prediction error in equ. (5. 17) tends only for limT P -* T(P) 
to zero. 

— The local deviation of the actual gravity field from the assumed iso- 
tropic and stationary global field causes an error AC in the covariances. The 
influence of AC on the prediction result will be of low order but considerable 
for the prediction error £ P (Rummel (1975)). 

It should be mentioned that the solution of Bjerhammar (1973) is not 
nonstationary. For nonstationary quantities related to a sphere it is impossible 
to find a covariance function in a closed form expanded into Legendre poly- 
nomials. 
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In the application of linear prediction methods we are faced with two problems : 

1. hi order to avoid filtering of desired frequencies the covariances have 
to be at least of the same degree as the output quantities. This means that for 
prediction of l°x 1° anomalies the degree of the covariances has to be higher 
than n=180. This causes serious trouble because gravity material necessary 
for a precise computation of covariances of such high degree is not available 
with sufficient density. The available terrestrial anomalies will lead to a not 
representative covariance function. The best possible covariance function 
based on the current gravity material is derived by Tscherning and Rapp (1974). 


2. For a high number of observations the autocovariance matrix C of 
equation (5.16) grows to a large dimension. The inversion of this large 
C -matrix is a serious problem comparable with the difficulties connected with 
the necessary inversion of the matrix in the deterministic approach. By 
applying stepwise methods the trouble can be reduced. 

5. 2 Deterministic approach 

The problems related to the nonstatistical solution for the downward con- 
tinuation are worked out in some detail by Schwarz (1971) for the planar approx- 
imation. The results derived there are valid with some modification for the 
spherical approach too and need not be repeated in detail. We have to 
solve the discrete version (5. 14) of an integral equation of the first kind (5. 1) 
for Ag P 7 . Restrictions of two types are usually not avoidable: 

The sample rate of the observations and their noise level allows only a 
limited resolution up to degree N. This restriction is not inconvenient because 
the frequency limitation is also necessary to stabilize the solution of the down- 
ward continuation problem. Because of the orthogonality relation between dif- 
ferent frequencies, the finite solution is consistent with the infinite one. The 
results up to degree N may be treated as reference for the evaluation from 
N+l to a higher frequency N+k, made possible by a more sensitive technical 
equipment. 

In addition f data will be usually gathered only in a limited area. When we 
assume the satellite moves in a near circular orbit then with a certain approx- 
imation the observations are given on an outer sphere with surface t . The sur- 
face of the sphere may be divided into a part T ( where observations are given 
and a part t 0 with no information, t= t, + ^ . On the other hand it is impossible 
to derive from the limited distribution of observations a global solution. 

A l a x 1° block global solution would require 64, 800 unknowns (= 44, 000 
equal area blocks). Thus also the surface a of the boundary sphere is divided 
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into an inner region CT| , where gravity anomalies will be derived aid a truncated 
outer region Ob , 0\ + a 0 = a. The part CT| will usually be somewhat larger than 
the subregion of t ( (the projection of t, onto <y). 

We rewrite equation (5. 14) as: 

(5. 18) T P? = Aj Ag Pj Agp Q » 

i 

where, for example A 0 transforms gravity anomalies from the domain Ob to 
disturbing potential values in the domain t, . The problem is always underdetermined 
because of the limited number of given T Pj - values. Only by smoothing, i.e. 
choosing the block size of the unknown gravity anomalies large the number of un- 
knowns decreases to a finite number m. If the number of T P) -values is higher 
than m we get from this smoothing procedure an "overdetermined" problem. 

For the finite solution with m unknown Ag-values in Oj , T P has to be 
corrected by AoAg (> 0 , equation (5. 18) becomes 

(5. 19) T P) - A 0 Agp Q = A' Ag P| 

When no gravity information on Ob is available we have to neglect the correction 
ApAgpo. This will lead to a systematic error in the results. The magnitude 
of this error depends on the size covered by o 0 , for lim Ob - 0 the solution con- 
verges to a global solution. 

If 0| is exactly the projection of T, onto a with respect to the common 
origin of both spheres the influence of A^ Ag P on the derived Ag P will grow from 
the center of t, to its border and result in a strong bias in Ag P _ . 

But it is also possible to build up every element of T P by a linear com- 
bination with Agp ( elements covering a certain cap, compare with Schwarz, 

(1971, p.40). Then oj will be larger than the projection of t, onto a and the 
error influence of ApAg Po will be about the same for every T P . Usually this 
type of solution is preferred. 

But we may also try to avoid or at least reduce the influence of Ao Ag P by 
computing it from available information in a 0 . Most precisely the outer part is 
taken into account by computing a correction A' 0 Ag P in Ob from all available 
terrestrial gravity anomalies and interpolated values in unsurveyed regions for 
all elements of Tp in satellite elevation. This procedure will lead to an optimal 
combination with gravity information especially when a t consists mainly of 
ocean regions. The large drawback is that the computation would be very time 
consuming. 

An approximation for the r.m. s. influence of A^Ag Pfl may be derived 
assuming a special mode for the degree variances v n (Ag) such as expressions 
shown by Rapp (1972) or by Tscherning and Rapp (1974). 
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Modifying the well known truncation formula of Molodenskii (Heiskanen 
and Moritz, p. 260) we get 


A I * * H 


/ 180 n 

/ St(r,0)Ag(p')da = 6T(P) 

jp= 0o a 


N 

r.m.s. 6T(P)=± y Q n 


V n+1 




V n (Ag) 


* 


( In any case the choice of a high reference field will cut down the influence 
of ApAgp o , for the remaining frequencies are more of local nature. For the 
solution of the discrete problem (5. 14) as well in the nonsingular as in the 
singular case, Bjerhammar (1973) gave some equations depending on the desired 
optimality criterion. The solutions are identical and equal to the prediction 
solition , equ. (5.7) for Ag and T assumed continuous on the spheres t and a. 

Bjerhammar's method for Ag^Ag^ = min. and A 0 with dimensions nxm, 
n<m (underdetermined) leads to: 

(5. 20) Agp ( = (A,' )Io T P , - (a|)Io AgP 0 

= (aJ)'[(A[) (Aj )'] _;L Tp, - (A| )'[a| )(AO'r 1 A^Agp o 

where usually the second term ( A_[ ) r [ ( a] )(A| /] -1 Ap Agp Q is neglected. The inverse 
(Aj ) 1 " 0 1 is the result of an iteration or approximation process, compare again 
Schwarz, (1971). 

Summarizing the problems for this type of solution we are faced with 

(1) A frequency limitation expressed by a certain block size which is caused by 
the limited resolution in the observation equipment and by the spacing between 
successive observations. 

(2) The error influence of the neglected outer zone Ob on the boundary sphere in 
the least squares solution. 

(3) The optimal choice of the inverse matrix (a[ )Jq . 

Now we try another solution for the downward continuation based on the 
inverse formula (5.4). We assume that there exists a downward continuation 
operator B which allows a continuation of the type 
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Ag(P 7 ) - B(I*, P)T(P)dT with P€T 

T 

In the global continuous limit such an operator exists in form of the inverse up- 
ward continuation operator developed into spherical harmonics, equ. (5.4). The 
discrete form of this operator shall build up the matrix B for the inverse problem 

(5.21) Ag,'=Bp' P Tp 

with 

N -> f n+l 

bp »V I [rJ (2n+l)P n (cos 0PjP^) 

for a certain element of B. 

Equation (5. 21) looks similar to the prediction solution, but in constrast 
to the prediction we try to approximate the operator directly from the continuous 
formula and not from the properties of the data expressed in covariances. This 
solution seems surprisingly simple, i.e. no matrix inversion seems necessary 
which with the same assumptions as before should hardly be true. Data Tp are 
only given in t, , and we get 

(5. 22) Agj> = Bj T| + Bo To 

where now Bo maps from the domain t 0 into the domain CT| . hi the solution we try 
to take into account the second term of equ. (5.22) BpT 0 . The elements of T 0 
for the zone t 0 may be expressed in terms of gravity anomalies on a by 

(5. 23) Tp = Aj Ag P; + Aq Agp 

Inserted into equ. (5.22) it will lead to 

( 5. 24) Ag Pi ' (I- B^ A° ) = b| Ti + BjJ A^ Agp Q 

and we have to solve again a difficult problem as before for Ag P) , as long 
as Bp and Aj are not orthogonal - which they never will become due to 
integral equation (5. 1), and as long as B<J or A? are not zero matrices -which 
is only true for the continuous and global limit where equ. (5.24) converges 
to equ. (5.14). 

By comparing these two methods expressed by equ. (5.20) and equ. (5.24) 
it becomes obvious that in both cases only the limitation of cr and not that of t 
has to be kept in mind. The limitation of T P to t, causes a high smoothing 
of the Solutions for the condition Agp Agp = min. 
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In order to make the problem in the prediction and in the least squares 
adjustment analysis more local a N l0 w degree reference field can be introduced. 
For the deterministic solution it reduces the influence of AT P = AoAg P „ in the 
least square prediction the residual covariances will decrease more rapidly to 
the zero line and then oscillate with (N lo w +l) zeroes, compare Meissl (1971), 
which diminishes the importance of the not available data of the outer zone. 
Because of the nonglobal distribution of all, up to now, derived satellite data 
and terrestrial measurements, the sets of spherical harmonic coefficients 
estimated from the limited information are disturbed by aliasing and dependence 
of the theoretically independent coefficients. Noise to signal ratios are rapidly 
increasing to 1:1 almost for low degree potential coefficients. The small con- 
fidence for the low degree field and the very systematic effect of aliasing 
caused by frequencies higher the Nyquist frequency, will result in a 
bias in the desired gravity anomalies based on this reference field. It is 
therefore necessary to accompany the local experiments by a global analysis 
to achieve a considerable improvement of the low degree and order harmonic 
coefficients. The simulation for the GRAVSAT/GEOPAUSE mission done 
by Koch-Argentiero (1974) indicates a realistic possibility for the low degree 
improvement. For the same experiment a local derivation of gravity anomalies 
may also be carried out by satellite to satellite tracking. 

The downward continuation problem was considered in thfc chapter in the 
form of equ. (5. 1) i.e. , as a computation of gravity anomalies on the surface a 
from the disturbing potential in satellite altitude, though for the discussed satel- 
lite techniques- satellite to satellite tracking and gradiometry-the gravity anom- 
alies shall be derived from the gradient VT of the disturbing potential or its ten- 
sor components. But this fact is not essential for the topic as a downward con- 
tinuation problem, considering the circumstance that the case formulated with equ. 
(5.1) is the worst possible. As shown in Table 2 and Chapter 4, the high fre- 
quency content of the quantities in consideration i.e. VT and gravity tensor, is 
more significant than in the disturbing potential itself. The problem will there- 
fore be more local. The operators connecting VT and the gravity tensor with Ag 
will not amplify uncertainties in higher frequencies to the same amount as the 
operator connecting T and Ag. The specific formulas necessary for the least 
square adjustment given the components of vt are explained in detail by Hajela 
(1974, p. 25)t those for the Ag computation from the tensor components of the 
disturbing potential by Reed (1973 p. 68). 

The application of the statistical prediction method for the analysis of 
satellite to satellite tracking should be divided into two steps. First, one may 
deduce the components of VT from satellite observations by least squares adjust- 
ment. Then from the gradient the gravity anomalies are computed by prediction. 
The partitbned analysis avoids the difficulties connected with the computation of 
auto-and cross-covariances for the sampled observations necessary in a direct 
one step prediction solution. 
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In addition, it is possible to estimate in the first step the VT - values only 
for a grid. This would reduce the matrices involved in the prediction step to a 
convenient size but on the other hand reduce the resolution. 

To predict gravity anomalies from VT means estimation of a scalar 
quantity from a vector process. Thereby we have to keep carefully in mind 
that orthogonal vector components for two distinct points P and Q are corre- 
lated even for a stationary and isotropic process. The covariance functions 
of a vector process becomes a tensor. 



Cu 

C l3 

Cl3 

C»(VT(P), VT(Q)) = 

C* 

^33 



Cax 

^33 

C 33 


It may be expressed in the Karman-Taylor decomposition, (Grafarend (1972)), 
as a function of the longitudinal and transversal covariances. The derivation 
of the tensor components for a vector process is explained by Moritz (1973, 
p. 60). The prediction formula (5. 15) is modified for the case of satellite to 
satellite tracking to 

(5.25) 


For the gradiometer experiment the prediction equations for the tensor com- 
ponents of the disturbing potential are derived in Moritz, (1970 p. 34) and 
need no further explanation. 

It should be mentioned that the prediction formulas are well suited for 
combination solutions with terrestrial and satellite observations of different 
type. For example the combination of satellite to satellite tracking data with 
the gradiometer data to take advantage of the specific wavelength properties 
of both methods does not have fundamental difficulties. 


6 . Proper Choice of the Inner Radius R 

In all considerations up to this point, reduction of the gravity data to a 
sphere a with radius R was assumed without specifying the semi-diameter in 
detail. Because the potential is regular and harmonic only outside the attract- 
ing masses also the continuation integral (3. 1) is valid only in the outer space. 

For the spherical formulation of the deterministic downward continu- 
ation problem a minimum semi-diameter for the boundary sphere fulfilling 
this condition would be R= 6384.403km due to Snowden (1968), (Rapp (1969)). 
The sphere would enclose all masses and touch the earth at a latitude of 
0= -1°28' in an elevation of 6272m. The advantage of reducing all gravity 
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quantities, i. e. satellite and terrestrial data to such an outer sphere- called 
geosphere by Bjerhammar (1967)- is mentioned by Eeg and Krarup (1973). 

They also explain the essential drawback: In a second step the derived poten- 
tial lias to be reduced from the geosphere to the actual surface of the earth 
where we finally need the information. Therefore we are again faced with all 
theoretical disadvantages of downward continuation as an "improperly posed 
problem". In addition, the reduction has to be carried out to a complicated 
surface oo"a closed smooth surface in £2 surrounding and arbitrarily near to 
the boundary of O" (Krarup, 1969) (£2. . . outer space of the earth). 

It complicates the solution essentially. More generalized one may 
interpret the problem as a free boundary value problem with additional physical 
and geometrical restrictions necessary for its solution, Grafarend (1972). But 
for its practical realization there is by far no ray of hope. 

The deviation 6 r for a point P E on the surface of the earth with geo- 
centric radius vector P to the geosphere with radius R - compare Figure 5 - 
reaches from zero to about 28km. 
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In satellite geodesy up to now only global solutions for a low degree field are 
computed. The introduced radius R is usually the 

equatorial radius a e = 6378. 16 km 

-compare for example Kaula (1969, equ, (5) and equ. (7)). Calculating with the 
derived coefficients gravity quantities by using formulas like equ. (2. 1) or equ. 
(2. 5) is nothing but a downward continuation of the gravity information to a sphere 
with radius a e . The drawbacks of this common practice are: 


1. A sphere with semidiameter a e does not enclose all masses. 
Therefore, for large areas the deduced quantities loose their physical meaning. 

2. The deviation 6 r of the earth's surface point P e to the point V* on 
the geosphere should be considered at least theoretically. The amount of 6 r 
for R = a e varies from -6km to +22km. For low harmonics the downward 
ccntinuation effect due to 6 r is not very high, as may be seen from Table 5a 
and 5b. The deviation in r. m. s. (Ag)^* caused by neglecting §R is for the 
gravity field up to ng = 20 about 1. 5% ana reaches 4% for some extreme areas. 
For the current accuracy level this influence may be neglected. For small 
block sizes the same deviation grows up to 10% and should no longer be ne- 
glected because of its strong systematic dependence on latitude and elevation. 




\ mgal 2 
altitude^\. 

var(Ag) 3>a0 

var(Ag) 3>3S var(Ag) 3>73 

- 5°x 5° ; =2.5 u x2.5° 

var(Ag) 3#ie0 
= l°x 1° 


235.7 

364.9 j 559.4 

1 

864. 6 

R - Rejo 
10 km 

227.8 

j 

346. 5 513. 1 

730 . 5 

R - Rej « 
22 km 

217.8 

323.8 459.5 

599.5 


Table 5a: Gravity Anomaly Variance Between Degree n x and n s 

For Low Altitudes 


varfAg)^ 


l M 

n, ^ J 


— mgal 3 ; 
n+B 


A = 425 . 28 mgal 2 , B = 24 
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n 

3-20 

3-36 



lkm to 10 km 

% 

n 

2.5 

4.2 

6.8 

lkm to 22km 

% 

3.8 

n 

9.3 

10.6 


Table 5b: Change of |r.m.s (Ag). | in percent from 

altitude 1km to altitude 10 or 22 km. 

This error source may be considerably reduced by introducing a 
spheroidal boundary surface instead of the geosphere as proposed by 
Hotine (1967). The computations for ellipsoidal harmonics are by far not as 
convenient but an ellipsoidal boundary surface will become necessary in the 
near future for several applications as for example high accurate geoid comp- 
utation, Lelgemann (1970). The transformation formulas from ellipsoidal to 
spherical harmonics are shown in Hotine, (1967, equ. (12) and (13)). 

A second and more promising approach is the computation of "gravity 
anomalies" at a boundary sphere completely embedded in the earth's surface 
with its origin in the earth's eenter-the so called Bjerhammar sphere. A 
possible choice for the Bjerhammar sphere would be 


R BJe = 6356.7km. (which is somewhat smaller than the semiminor axis 
of the earth) 

The application of this method is today generally recommended, Hotine 
(1969), Krarup (1969), Moritz (1970). Although the physical meaning of the 
reduced values is lost it is in best agreement with the practical restriction of 
an always finite number of measurements, to a second step the deduced 
'gravity anomalies' have to be upward continued and can be compared with 
the terrestrial gravity information. 

By applying the least squares prediction approach the question of a 
proper choice of R seems not to be very troublesome. The geocentric radius p 
to a surface point P e may be inserted into the necessary covariance C-Ag P t ? 
expressed by equ. (5.9a) to lead with equ. (5. 7) immediately to the desired 
boundary value. It is a very convenient procedure, but it does not ensure con- 
vergency as long as the convergency of the applied covariance function on the 
surface of the earth is not proved. When in the covariance expressions the 
empirical degree variances are replaced by practical applicable analytic ex- 
pression a fixed radius has to be found. A detailed analysis for a proper 
formula was done by Tscherning and Rapp (1974). The introduced expression 
is 
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(6.1) var(Agp') = V s n+ — 

' ’ ' SP ; ^ (n-2)(n+B) 


n=3 


with the unknown parameters s, A and B where 


s = 


Rl 

Rp/ 


R . . . radius of an inner reference sphere (Bjerhammar sphere) 

R P ' ... central distance for a surface point p', for the analysis to be replaced 

by 

R e . . . mean radius of the earth 


The result of the least square fit is (Ibid, p. 22, Table 7) 
s = 0.999617, A = 425. 28 mgal 3 , B = 24 and gives R = 6369. 8km. 

A Bjerhammar sphere with this semi-diameter would not be embedded 
in the earth and makes another choice for R necessary for downward continu- 
ation applications. For our problem the fitting procedure shall lead to a regular 
and harmonic expression for the covariance function that converges on the 
Bjerhammar sphere with minimum variance and least squares deviation to the 
empirical data at the earths surface. The theoretical formulation of such an 
adjustment for the potential was done by Krarup (1969, p. 54) based on the 
Runge's theorem. The analysis for equ. (6. 1) based on the semi diameter R BJe = 
6356. 7km gives the results, using the same procedure and data as in Tscherning 
and Rapp (1974). 

s = 0 . 995516, A= 705. 094 mgal 2 , B = 43. 

The covariance expression for downward continuation based on these values is 


< 6 ' 2) (^fj 


The drawback of the result lies in the fact that the variance of the point anomalies 
is too low (var(Ag) 3t500 = 1081. 7 mgal 3 ) and will therefore lead to an additional 
smoothing in the least squares prediction solution. Especially for high elevations 
the covariance function decreases very rapidly with increasing degree. In 
addition for all expressions of the form of equ. (6. 1) convergency on the 
Bjerhammar sphere cannot be ensured. Unrealistic amplifying of the high fre- 
quency content in the predicted quantity near or on the Bjerhammar sphere may 
be the result. 
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tli th. 

The quotient Q of the (n+1) term a n+1 and of the n m term a n of 
var(Agp'), equ. (6.1), withs=l leads to 

a n+ i _ (n a - 2n) (n+B) 

^ a „ (rr + 2n+l) (n+B+1) 

and has the limit 

lim Q = 1, 

n -* “ 

which should be smaller one for convergency (d'Alembert criterion). Therefore 
expanded analysis with more complicated analytical expressions for var(Ag) ful- 
filling the demanded properties would be very usefiil. 

7. Conclusions 

Any derivation of surface gravity information from satellite observations 
is a downward continuation procedure. If the downward continuation is carried 
out in a global solution by means of a spherical harmonic expansion of the de- 
sired gravity quantity to a sphere containing the masses of the earth (geosphere) 
no essential difficulty arise but the values related to the geosphere are because 
of their distance to the earth's surface or the geo id, of no use for practical appli- 
cations. 

Therefore, in practice the spherical harmonic coefficients are usually 
derived in a satellite method with respect to a sphere with radius a e (a e ... semi- 
major axis of the earth). A sphere with radius a e intersects the earth's surface 
in some areas, and in most areas lies outside of it. Thus, the disturbing potential 
T derived from this set of spherical harmonic coefficients by 


N n 

T(P)=Y j £ (c».(T)R„(P') + s na (T)S na (P')), P'€CT 

n = 2 »= 0 

( a yn+1 

~| ■ = 1 , or r= a e 

has in a strict sense no physical meaning because the solution of the boundary 
value problem is valid only for functions harmonic and regular outside the boun- 
dary surface. The resulting problem of divergency for the spherical harmonic 
expansion at and inside the surface of the earth becomes not obvious in a low 
degree solution, but the physical meaning of gravity quantities calculated from 
this set of spherical harmonic coefficients is doubtful. These considerations 
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get practical relevancy in a combination of satellite derived gravity values with 
terrestrial gravity data, Rapp (1969). 

Because of a strictly analytical continuation of gravity quantities to the 
actual surface and its interior is not possible, (Moritz (1961), Bjerhammar (1969)), 
the purpose of physical geodesy may be reformulated to the "determination of a 
gravity field that is compatible with the given discrete observations", (Heiskanen 
and Moritz, p. 321). This formulation is in agreement with the real situation 
where always only a finite number of observations is given. 

W' tb. this more realistic statement we may look at the local derivation of 
gravity anomalies from satellite to satellite tracking or satellite gradiometry 
rather optimistically. For this approach we are theoretically faced not only 
with the question of a proper choice of the boundary surface but also with the 
restriction of the boundary surface to a finite region where we try to derive 
gravity information. Results of simulation studies of Reed (1973) for satellite 
gradiometry and Hajela (1974) for SST seem to be promising. But as any sim- 
ulation they suffer from some simplications which shall be tept in mind for 
planned practical experiments: 

a) For simulation studies the input gravity anomalies are thought to be referred 
to a sphere with earth's mean radius 6371km. The gravity anomalies resulting 
from the simulation are again related to the same sphere and can easily be 
compared with the input data. In a practical experiment the derived gravity quan- 
tities related to the boundary sphere have to be transformed to be comparable 
with terrestrial gravity material. 

b) The simulation data are generated by mean gravity anomalies of a certain 
block size. White noise is added to simulate observation errors. In reality the 
satellite is affected in its motion by the whole gravity spectrum. The obser- 
vations are done with a certain amouti/t of normal distributed white noise and 
with some systematic errors. The error level defines the limit of the resolution 
in the frequency range. But this limit is very flexible. Therefore in a certain 
frequency range the noise will be somewhat biased and cannot be filtered out by 
least squares methods. In practical experiments, this bias will affect the re- 
solved spectral content of the derived mean gravity anomalies. 

c) The downward continuation problem expressed by integral equation (5. 1) is in 
any way undertermined, because we try to derive a continuous function on the boun- 
dary surface from a finite number of observations. 

For practical applications we approximate the integral equation by a system 
of linear equations with m unknowns for n observations. The three different cases 
m<n, m = n, and m>n express thereby the degree of approximation of the integral 
equation by the set of linear equations. For a practical situation a solution with a 
certain m and n, and m<n may be optimal in terms of optimal mean anomaly re- 
covery. From the theoretical point of view, with equ. (5. 1) in mind, it does not 
express an overdetermined least squares adjustment problem. 
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d) The local solution is in contrast to integral equation (5. 1) or its finite 
approximation by linear equations which assume global recovery of anomalies. 

In the simulation studies the satellite "observations" are generated from gravity 
anomalies of a finite region at the surface of the earth whereas in reality a 
satellite is influenced in its motion by the global gravity field of the earth. 

Gravity anomalies are then derived in a somewhat smaller region from the 
simulated observations. Therefore a numerical analysis of the influence of the 
outer zone on the derived gravity anomalies would be of great use. 

Based on these considerations we propose the following solution for the 
local recovery of gravity anomalies which is in best possible agreement with 
theoretical demands and practical restrictions. We may choose between two 
different approaches , one is the least squares adjustment or inversion procedure 
based on the approximation of the integral equations for the VT - components and 
the disturbing potential tensor-components^ the other is the least squares prediction 
method. For both types two steps are necessary in their solution. 

Inleast squares adjustment we first derive "gravity anomalies" with no 
real physical meaning at a Bjerhammar sphere totally embedded in the surface 
of the earth by downward continuation.by such procedures as describedby Reed (1973) 
or Hajela (1974). Then by upward continuation these quantities are transformed to the 
surface of the earth, where the resulting anomalies willbe compared with terrestrial 
gravity anomalies. The second step, the upward continuation procedure is not connected 
with principal problems. 

In the least squares prediction solution we derive first from SST -obser- 
vations or gradiometer observations gravity quantities in satellite altitude- 
the components of VT r in SST and the independent tensor components Tx , 
i= 1, 2, 3 and j = 1, 2, 3 in gradiometry. This intermediate step is necessary be- 
cause we have no information about the covariance functions for the observations 
and about the crosscovariances connecting the observations with the gravity 
quantities. In a second step we derive with downward continuation by least 
sqiares prediction immediately point-or mean anomalies at the surface of the 
earth, ( Moritz (1970), (1973)), that may be compared with terrestrial anomalies. 
Practical problems related with both procedures are described in Chapter 5. 
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